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On the basis of the Keldysh method of non-equihbrium systems, we develop a theory of electron 
tunneling in normal-metal/superconductor junctions. By using the tunneling Hamiltonian model 
(being appropriate for the tight-binding systems), the tunneling current can be exactly obtained in 
terms of the equilibrium Green functions of the normal metal and the superconductor. We calculate 
the conductance of various junctions. The discrepancy between the present treatment and the well- 
known scheme by Blonder, Tinkham, and Klapwijk is found for some junctions of low interfacial 
potential barrier. 
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I. INTRODUCTION 

One of the powerful methods detecting the quasipar- 
ticle states in a superconductor is to measure the con- 
ductance of a junction made up of a normal-metal and a 
superconductor (NS). There have been developed many 
theories describing the electron-tunneling phenomenon. 
In the case of the high interfacial potential-barrier limit, 
the linear-response theory is a well-known description.^ 
But it is not valid for describing the electron transport 
in the low potential-barrier limit. 

To calculate the conductance in more general cases. 
Blonder, Tinkham, and Klapwijk (BTK) have developed 
a theory by supposing that the system is in such a non- 
equilibrium state that only the incoming particles have 
equilibrium distributions.^ This theory has been widely 
used for analyzing the tunneling phenomena in various 
NS junctions, and also has been extended for investigat- 
ing electronic tunneling in Josephson junctions.'^ 

When a finite voltage is applied to a junction, the 
electron transport in the junction is a non-equilibrium 
process. We would like to consider the case when the 
current passing through the junction is a constant. The 
electron transport process is then a steady state. Such 
a non-equilibrium problem can be solved by the Keldysh 
approach.^ In fact, this approach has been applied by 
a number of investigators for studying the tunneling in 
junctions of normal metals and the electron transport 
under impurity scattering.'' 

In this paper, we present a tunneling theory along this 
direction. We will start with a tunneling-Hamiltonian 
model defined in a square lattice. This model is appropri- 
ate for the tight-binding systems. The tunneling current 
can be exactly obtained in terms of the equilibrium Green 
functions of the normal metal and the superconductor. 
By so doing, all the effects of external voltage on the 
tunneling current can be rigorously taken into account. 
Moreover, it can be extended to study the tunneling in 
the point-contact junctions as in the scanning-tunneling 
microscope measurement. 



II. FORMALISM 

We consider a junction consisting of a normal metal 
on the left side and a superconductor (SC) on the right 
side. In the Nambu representation, the tunneling Hamil- 
tonian describing the electron-transport processes in the 
junction is given by 

Ht = ^^{clfrlCl + cjfirCr) (1) 
Ir 

where cj, — {clp Cri) is the field operator for particles in 

the right superconductor, and cj is similarly defined for 

the left metal, fri = fl = iodj/r - yi\)(^3, and y,. and yi 
are respectively the coordinates of the sites r and I along 
the interface. The r and I summations in Eq.(l) run over 
the edge (interface) sites on the two sides of the junction, 
respectively. The function to{\yr — yi\) may be taken as 
real. For simplicity of description, we suppose that the 
lattice sites {r} along the edge are equally spaced as the 
same as {I}. Suppose there is a voltage V apphed be- 
tween the junction, the total Hamiltonian of the system 
is given by 

H = + Ht ^ Hi - eVNi + Hr + Ht, (2) 

where Hi and Hr are the intrinsic Hamiltonians of the 
left metal and the right superconductor, respectively, and 
Ni is the total electron number of the left metal. We here 
adopt the tight-binding model for Hr which contains a 
hopping term and an attraction term. For Hi, we keep 
only the hopping term. 

To define the tunneling-current operator, we first con- 
sider the charge operator for the right SC. Apart from a 
constant, it can be written as 

Q^-e^clasCr. (3) 

r 

The operator of current through the junction from left to 
right is then obtained as 

i ^i[H,Q] = ie ^(cJ.craTwQ - clfirCrsCr). (4) 

Ir 
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Now, let us choose the unperturbed state described by 
Hq as our reference system. This reference system con- 
sists of the unperturbed normal metal and the SC on 
two side of the junction, each of them in its own equi- 
librium state. For the purpose of employing the grand 
canonical ensembles, we use Ki — Hi — [jii + eV)Ni and 
Kr = Hr — firNr to describc the normal metal and the Sc, 
respectively. Here, and Hr are respectively the chem- 
ical potentials of the normal metal and the SC, and Nr 
is the total number of electrons in the SC. At the steady 
state, we have Hr = Hi+ eV in order to maintain charge 
neutrality in the bulk of each side. To calculate the sta- 
tistical average of a physical quantity, we need to write 
the related operators in the interaction picture. An op- 
erator of physical quantity, e.g., the current I{t), in the 
interaction picture at time t is defined as, 

i{t) = exp(ii/°t)/exp(-ii74). 

This operator can be further rewritten in terms of the 
field operators, 

i{t) = -2eImY^ 4{t)a3fri{t)ci{t), (5) 

Ir 

where cl{t) = exp{iKrt)cl exp{—iKrt) (and a similar def- 
inition for ci{t)), fri{t) = flit) = friexpiieVtaa). The 
form for I{t) as given by Eq. (5) is convenient for the 
statistical average over the grand canonical ensembles. 
Similarly, the tunneling Hamiltonian can be written as 

Hrit) = J2[4{t)Trl{t)ci{t) + cl{t)fir{t)Cr{t)]. (6) 
Ir 

For applying the Keldysh method, it is convenient to 
define the field operator, 

4it) = [cl{U),cl{t.)] (7) 

where the subscripts -|- and - on time t means the 

operators defined in the time branches (—00,00) and 
(00,— 00), respectively. Accordingly, we define a pertur- 
bation Hamiltonian, 

Hc{t) = Y,[4im'imit) + 4mut)Mt)] («) 

Ir 

where 

The matrix is the third Pauli matrix defined in the 

space corresponding to the two time branches. To distin- 
guish with that, we reserve 173 as the third Pauli matrix 
defined in the particle-hole space. The Green function is 
defined as 

Gij{t,t') = -i{T[S,Mm]{t')]) 



/oo 
dtHc{t)] 
-00 

where T is the Keldysh time-ordering operator. 

With the above definitions, the current under the sta- 
tistical average can be expressed as 

I = e^Re'TTasfri{t)Gir{t,t), (10) 

Ir 

To calciilate the ciirrent, we need to know the Green 
function Gir{t,t). It can be determined from the Dyson 
equations. 

Let L and R denote the Green functions (as 4x4 ma- 
trices) for the left metal and the right SC, respectively 
(with the superscript for the unperturbed ones). By 
assuming that the system is uniform along the direction 
parallel to the interface, we can then work in the mo- 
mentum space. Here, the momentum is parallel to the 
interface. The Dyson equations are 

Gk{t,t')= [ dhLl{t,h)T^\h)Rk{h,t/) (11) 



Rk{t,t') = Rl{t,t')+ / dti / dt2RUt,h)^k{ti,t2)Rk{t2,t') 



Mtut2) = mh)Ll{tut2K\t2) (13) 

where T^{t) = TzTk expiieVtas) , Tk = to{k)a3, and the 
range of time integrals is from —00 to 00. Note that 
the Green function L^{ti,t2) = L^{ti — ^2) consists of 
four diagonal matrices. The factors exp(ieyiiCT3) and 
ex.-p(ieVt203) commute with the matrix Lff.{ti,t2)- The 
self energy Sfc(ti , t2) = J^k{ti—t2), and thereby the Green 
function Rk{t, t') = Rk{t — t') are functions of time differ- 
ence. We can therefore take the Fourier transformation 
of the Dyson equations. In the frequency space, these 
equations have the usual forms except 

S,(a;) = nmliu^ + eVas)T^,\0). (14) 

With the help of the Dyson equations, we can write the 
factor Tri{t)Gir{t,t) in the expression of / as 

/oo 
—T,i:k{io)Rk{iv). (15) 

Inserting Eq. (15) into Eq. (10) and taking the trace of 
time-branch space, we have 

/oo J 
— i^ReTr<T3M+(L/ii° + L\Rf)M_ 

K 

(16) 

with 

M± = [l-tlLlRl]-\ 
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Lf= tanh[(a; + eFc73)/2fcBT](L° - L°_), 
Rf = ta.nh{u)/2kBT){Rl - R°_), 
= i^ot = u> + eVaa + iO), 

Rl = R°J = R°{k,w + iO). 

Here and (L^ and i?" ) are the retarded (ad- 
vanced) Green functions (as 2 x 2 matrices in the Nambu 
space) of equiUbrium state, Lf and Rf are the Keldysh 
functions, = |io(fc)P, ks is the Bohzmann con- 
stant, and T is the temperature of the system. By 
noting the relationships _R+(fc, — w) = —U2R-(k,ijj)u2, 
Lj^(k,—uj-\-eVaz) = —(J2L-{k,LO + eVaz)'^2, it is enough 
to only take the frequency integral in Eq. (16) in the 
range (0, oo). The front factors in the Keldysh functions 
take part of the roles of quasiparticle distribution func- 
tions. The additional term —eVa^ reflects the chemical 
potential shifts of the quasiparticles in the left metal. 



III. GREEN'S FUNCTIONS OF THE 
EQUILIBRIUM STATE 

To calculate the tunneling current /, we need to know 
the Green functions and If we know the wave 
functions tpn and energies En of the quasiparticles, e.g., 
for the SC, we can obtain i?° by 

R\k,.) = Y.^^ (17) 

n " 

where ipn takes the edge value. Since we have taken the 
Fourier transformation for the dependence on the coor- 
dinates parallel to the interface, the wave function ipnij) 
depends on the x-coordinates (normal to the edge) of 
lattice sites, j = {1, 2, • • •}; the edge value is tpni^)- 

For illustration, we here consider a d-wavc SC and sup- 
pose that the order parameter is constant everywhere. 
The wave functions can be determined analytically by 
the BdG equation. As an example, we consider the tight- 
binding model defined in a semi-infinite square lattice 
with a {11} edge. The BdG equation reads ^ 

J2HMj)=Eni^n(t), (18) 

j 

where Hjj = — /x(J3, -ffjj-i = —2tcoska3 — i2Asmkai 
for j > 2, Hjj+i = — 2fcosfcf73 + i2Asinfccri, otherwise 
Hij = 0, t is the hopping energy of electrons between 
nearest-neighbor sites, and A is the order parameter. 
Here, we have used the unit v^/o (with a the lattice 
constant) for the momentum k, and k is confined to a 
Brillouin zone (— 7r/2, tt/2). There are two kinds of solu- 
tions to Eq. (18): The continuum states and the surface 
bound states. 



The contimmm states are generally degenerate. To 
distinguish them, we can consider each eigen wave func- 
tion contains a unique incoming wave component or a 
unique outgoing wave component. We then characterize 
the wave function by the incoming wave number q^i or 
the outgoing wave number q^. For example, the wave 
function and energy of state can be written as 

Mi) = i^kAi) - E «^aV'L(i)]/^^' (19) 

a 

Ek^^ = ±E{q^,k) = ±^e2(g^,fc) + A2(g^,fc), (20) 

where ip^'s are the plane- wave solution to the infinite sys- 
tem, e{q, k) = —Atcosqcosk — fi (with /z the chemical po- 
tential), A((7, fe) = — 4A sing sin fc. The coefficients a^a 
are determined by the boundary condition at j = 1. The 
summation over a in ccj. (19) runs over all the outgoing 
components with E{qa,k) = E{qfj_,k). It is worth notic- 
ing that sometimes we may have complex q^ 's, the sum- 
mation then should be taken at those q^s corresponding 
to decaying waves. 

The number of the bound states is determined by the 
Levinson theorem.^ Under the assumption that the order 
parameter is constant, we only have the state with En = 
for each |fc| < km {km is very close to the Fermi wave 
number). For £"„ = 0, it can be shown that the two 
components Uk{j) and Vk{j) satisfy the relation, 

Vk{j) = iXuk{j), A = ±l. (21) 

Suppose Uk{j) = with z (|z| < 1) a complex quantity 
for the general solution. Corresponding to z, we have a 
complex number q = —i \og{z). The equation E{q, k) = 
determining the eigenvalue reduces to 

t{z + z~'^) cos k + X{z- 0"^) A sin k + iJ,/2 = 0. (22) 

The solutions to Eq. (22) are 

z± = [-M ± - (cf - ci)]/(ci + Ac2), (23) 

where ci = 4t cosfc and C2 = 4Asin/c. Note z+z_ = (ci — 
Xc2)/{ci+Xc2), therefore A = sgn(fc) whereby \z^Z- \ < 1. 
The wave function is given by 

Ukij) = (4 - zi)/Nk, (24) 

with Nl = 2[(1 - |^+|2)-i + (1 - |2_|2)-i - 2Re(l - 
0^2-)"^] the normalization constant. This wave function 
satisfies the boundary conditions at j = 1 and j — > oo 
provided |z±| < 1. If /z^ < (cf — c|), then z+ and Z- 
arc complex conjugates of each other, and \z±\ < 1. On 
the other hand, if /x^ > (c^ — c|), both of them are real. 
In this case, there may be no bound state unless both 
\z±\<l. 

With the knowledge of the wave functions, the Green 
function i?° can be calculated by Eq. (17). As for L° of 
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the normal metal, it contains only the continuum states. 
The wave functions can be obtained immediately from 
Eq. (18) by setting A = 0. The resulted Green function 
is given by 

L\kM = - r dq^^p—. (25) 

IV. COMPARISON WITH THE BTK THEORY 

Obviously, the present treatment is a non-pcrturbative 
theory. It takes into account all the effects of the volt- 
age within the model. At this point, it is instructive to 
compare our theory with the BTK theory. In the BTK 
model, only the incoming particles in each side of the 
junction are described by the equilibrium distributions 
with the chemical potential shift of the left metal due to 
the external voltage. But, the outgoing particles are not 
described by the equilibrium distributions. The quasipar- 
ticle states in the whole system are determined by the 
Bogoliubov-de Gennes equation that is independent of 
the external voltage. The tunneling current is calculated 
as the result of the current by the incident particles from 
the left metal minus that from the right SC. In contrast, 
by the present consideration, the particle distributions 
are referred to the reference system. Since in the interac- 
tion picture, the tunneling Hamiltonian depends on time, 
there cannot be quasiparticle states for the whole system. 
Each state in both sides of the junction has its lifetime 
because of the non-equilibrium process between the inter- 
face. From the Green function, the lifetime of a quasipar- 
ticle is determined by the inverse of the imaginary part of 
the sclf-encrgy. In this approach, the transport process 
is treated by the equivalent of time-dependent perturba- 
tion theory to all orders, which leads to lifetimes. The 
electron transport is the process of quasiparticles decay- 
ing. On the other hand, in the BTK model, the transport 
process is treated by the time-independent perturbation 
theory to all orders, which determines the quasiparticle 
states in the whole system, with infinitive lifetimes for the 
continuum states. Therefore, the mechanisms of electron 
transport through the junction by the two theories are 
very different. 

For numerical comparison, we need to present the BTK 
scheme in the lattice model. The basic work in the 
scheme is to solve the BdG equation for the wave func- 
tions of quasiparticles in the whole system. An cigen 
wave function characterized by an incoming wave in the 
left metal can be written as the incoming wave plus all 
the reflected waves (including the Andreev and the ordi- 
nary reflections), with the transmitted waves in the right 
SC including all the outgoing waves. One needs only then 
consider the boundary condition at the interface barrier. 
By denoting the wave functions in the left and right sides 
respectively by ipiii) with j = { — 1, —2, • • •} and ipr{j) 
with j = {1,2,---}, the BdG equation at the interface 



barrier reads 

if_l,_2VK-2) + if-l,-lVi(-l) + T'fcVr(l) = Ei>i{-\), 

(26a) 

TfcVK-l) + -H'i,iV'r(l) + ifi,2V'r(2) = Eil,r{l)- (266) 

Eqs. (26a, b) are nothing but the boundary conditions. 
With the wave functions, one can immediately calculate 
the tunneling current according to the BTK theory. 

To see the difference between the present and the BTK 
theories, we have carried out the numerical calculations 
of the tunneling conductance 



for normal- metal/d- wave superconductor junctions with 
{110} and {100} interface at various barrier strengths. 
For presentation, we normalize G by Nc'^/tt {Ti = \) 
with N the total number of the lattice sites on one 
side of the interface. The basic parameters for the 
SC are, t = 176meV, hole concentration 6 = 0.15, 
attractive potential between the nearest-neighbor sites 
V = VlAmeV . The transition temperature Tc and the 
order parameter Aq are obtained as Tc ~ QQK, and 
Ao = 4A|t=o = IQ-lmeV, respectively. As being stated 
before, the Hamiltonian of the left metal contains only 
the hopping term. We assume that the hopping energies 
of both sides of the junction are the same. For simplicity, 
we choose tunneling matrix element as to{k) = to. 

The numerical result for the normalized conductance 
as function of V for an NS (d-wave) junction with {100} 
interface at T = is shown in Fig. 1. The tunnel- 
ing parameter to/t = 0.5 is used. Though the interfa- 
cial potential barrier at this parameter is not very high, 
the agreement between BTK and the present theories is 
very good. A small to means a high interfacial potential 
barrier. At the high potential barrier limit, both theo- 
ries reproduce the linear response result [1] . However, at 
to/t = 1 corresponding to a weak barrier, the discrepancy 
is clear as shown in Fig. 2. At weak barrier and small 
voltage \eV\ < Aq, the Andreev reflection is the predomi- 
nant contribution to the conductance in the BTK theory. 
Under the present assumption, however, the transport is 
due to the decay of quasiparticles in both sides. Such a 
decaying process is more complex than the BTK picture. 
The difference between the two theories at small \eV\ is 
mainly due to the different treatment of the tunneling 
Hamiltonian (i.e., time-dependent vs time- independent 
perturbation theory). The voltage effect in is impor- 
tant only at large \eV\, because the relevant dimension- 
less parameter is the ratio eV/Ep (with Ep the Fermi 
energy of the left metal). The voltage effect is more ev- 
ident at y < than at y > 0, because more precisely 
the parameter is actually \eV/{Ep + eV)\. At negative 
voltage, the chemical potential of the left metal shifts up- 
ward, resulting in electrons right below the Fermi surface 
within the energy range {Ep + eV, Ep) transferring into 
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the right SC. At positive voltage, the states in the energy 
range (Ep,Ep + eV) in the left metal are available for 
the electrons in the right SC to transfer in. 

In Fig. 3, we show the results for the junctions with 
{110} interface Sit to/t = 1. In this case, the results by 
both theories are in excellent agreement. The agreement 
is even better at smaller to. At \eV\ < Aq, the conduc- 
tance G is given by a broadened zero-bias peak. Actually, 
there are zero-energy bound states in the right SC near 
the interface, with lifetime due to tunneling. The tunnel- 
ing current is predominantly conducted by these states. 
The width of the broadening is mainly determined by 
the tunneling parameter to rather than by the external 
voltage. Because of the existence of these states, the par- 
ticle transmission through the junction for \eV\ < Aq is 
a resonant process. These resonance states exist in the 
BTK model as well. At least at eV — 0, both theories 
produce the same resonance states with the same energy 
broadening. Therefore, we can understand the excellent 
agreement near eV = 0. 

The discrepancy between the two theories is even more 
clear for the normal- metal/conventional-superconductor 
jimctions. Fig. 4 shows the result for an NS (s-wave) 
junction with {100} interface at T = and to/t = 1. 
The parameters for the SC are, the chemical potential 
H = —0.3t, the on-site pairing parameter Aq = 0.02t. 
The conductance predicted by the present theory is only 
about 78% of that of BTK for |ey|/Ao < 1 where the 
conductance is almost a constant. Qualitatively, the elec- 
tron transport in this junction is similar as that in the NS 
(d-wave) junction with {100} interface. The explanation 
for Fig. 2 applies here. 

V. AN APPROXIMATION SCHEME 

When eV/Ep <^ 1, the dependence of on the ex- 
ternal voltage is very weak. We can then drop eV in L*^. 
By this approximation, the conductance G is given by 

/OO 7 
^tlTvlm{R°_M_a3M+)a3lmL+g 
-°° 

(28) 

where g = eosh ^[(w + a3eV)/2kBT]/2kBT is only fact 
which depends on eV. In Fig. 4. the result by Eq. (28) 
is also plotted. At small voltage, the approximation is 
in very good agreement with our main theory. However, 
at large voltage, the approximation reproduces the BTK 
result. This clearly shows that the discrepancy between 
our main theory and the BTK theory at small voltage 
is not due to the voltage effect in In the case of 
eV/Ep <C 1, Eq. (28) is a simple but good scheme for 
calculation of the conductance. 



VI. SUMMARY 

In summary, on the basis of the Keldysh approach, we 
have developed a theory of electron transport in normal 
metal/superconductor junctions to all orders in the ap- 
plied voltage and the barrier strength. In the present 
scheme, the tunneling current is given in terms of renor- 
malizcd Green functions of a steady state. It can give a 
reliable description of the electron tunneling, including 
the ballistic transport in NS junctions. We have calcu- 
lated the tunneling conductance for various NS junctions 
using the present formalism and have compared it with 
the BTK theory. In most cases, both theories agree with 
each other. However, for some junctions of low barrier 
strength, the discrepancy between the two theories can 
be sizable. 
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FIG. 1. Conductance G as a function of the Voltage V for 
an NS (d-wave) junction with {100} interface at T = and 
to/t = 0.5. The present calculation (solid line) is compared 
with the BTK result (dashed line). 



FIG. 3. Conductance G as a function of the Voltage V for 
an NS (d-wave) junction with {110} interface at T = and 
to/t — 0.5. The present calculation (circles) is compared with 
the BTK result (squares). 
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FIG. 2. The same as Fig.l but at to/t = 1. 



FIG. 4. Conductance G as a function of the Voltage V 

for an NS (conventional SC) junction with {100} interface at 
T = and to/t = 1. The present calculation (dotted line 
with circles) is compared with the BTK (dashed line), and 
the approximated (solid line) results. 
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